Electrons in Dry DNA from Density Functional Calculations 
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The electronic structure of an infinite poly-guanine - poly-cytosine DNA molecule in its dry A- 
helix structure is studied by means of density-functional calculations. An extensive study of 30 
nucleic base pairs is performed to validate the method. The electronic energy bands of DNA close 
to the Fermi level are then analyzed in order to clarify the electron transport properties in this 
particularly simple DNA realization, probably the best suited candidate for conduction. The energy 
scale found for the relevant band widths, as compared with the energy fluctuations of vibrational or 
genetic-sequence origin, makes highly implausible the coherent transport of electrons in this system. 
The possibility of diffusive transport with sub-nanometer mean free paths is, however, still open. 
Information for model Hamiltonians for conduction is provided. 



INTRODUCTION 

The idea of using DNA chains as molecular wires in 
nano-devices has been proposed and is being explored in 
the nano-science community [Q, |^, || , given the extremely 
high selectivity of a DNA strand to pair with its comple- 
mentary partner, which facilitates the design of complex 
circuitry at the nanometer scale. Even though some re- 
searchers have chosen to coat |^ or modify Q the DNA 
wires in order to improve its conduction properties, the 
original idea was to use pure DNA, assuming that elec- 
tric current could flow along the wires. This assumption 
has proven to be less trivial than expected, however, and 
has originated an important controversy. 

Direct measurements show from low resistance DC 
conduction 01 (even superconduction |^), to no conduc- 
tion at all [p, \ including results with non-ohmic char- 
acteristics [pl |l(| or appreciable ac conductance . The 
transport has been proposed to be band-structure-like , 
by tunneling |12 , solitonic fl^, polaronic [|l0[ or 
facilitated by fluctuations |[lq , [l7|| . This rich scenario 
is not surprising when considering the comparably small 
amount of clean-cut information available (mostly the- 
oretical) about basic aspects of the electronic structure 
of dry DNA ||, (in all the experiments 

mentioned above the samples were dried in one way or an- 
other). In addition to the fundamental complexity of the 
system (large system size, weak interactions, structure 
floppiness), there is a substantial lack of detailed knowl- 
edge about the environment the DNA molecules are left 
with after the drying process. Most likely, residual water 
molecules will stay adhered to the chain, together with 
the cations counterbalancing the negative charge of the 
molecule, all quite uncontrolled. 



A theoretical simulation benefits from absolute control 
at the atomic scale. In this particular case, the system 
can be prepared in the best possible way for conduction. 
Previous simulations |2|, || have offered very 

insightful information for related problems, in spite of 
being based on relatively drastic approximations in the 
electronic structure, due to the system size. They were 
mostly aimed at the description of DNA properties in live 
conditions, as were the recent first-principles studies in 
wet DNA 0, |2i, I25II. 



Recent developments in linear-scaling first-principles 
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|27| , |28| methods based on density-functional theory 
[2g (DFT) allow the study of an infinite DNA double 
helix from first principles. DFT functionals based on the 
Generalized Gradient Approximation (GGA) have shown 
to provide very satisfactory accuracy in molecules akin to 
DNA |30, 31|. In this work, an extensive study is pre- 
sented on a set of 30 nucleic base pairs that validates 
the theoretical method. Then, an infinite DNA chain 
of trivial sequence is studied: all guanines (G) on one 
strand, all cytosines (C) on the other, in the Watson- 
Crick pairing arrangement. We simulate it completely 
dry and in vacuo, and in an A-helix structure, the most 
relevant one for these conditions. The electronic struc- 
ture is analyzed after relaxing the structure with the first- 
principles forces. Genetic sequence change effects are ex- 
plored on a similar system where one guanine-cytosine 
pair is swapped every eleven base pairs, then relaxed and 
analyzed. Preliminary results of this study were pub- 
lished in a joint experiment-theory paper Q. 
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METHOD 



STUDY ON BASE PAIRS 



The first-principles DFT calculations reported here 
were performed using the Siesta method |2^, 
which has already been applied to a large variety of sys- 
tems [Q, including biomolecules |3^. It is a fully 
linear-scaling method, both in the building of the DFT 
Hamiltonian and in its solution. GGA |^ was used 
for electron exchange and correlation. Core electrons 
were replaced by norm-conserving pseudopotentials | |3^ 
in their fully non-local formulation . A uniform mesh 
with a planewave cutoff of 150 Ry was used for integra- 
tions in real space |2^. The basis sets used are made 
of numerical atomic orbitals of finite range and are de- 
scribed below. 



The unit cell of the A-helix studied here includes eleven 
poly(G)-poly(C) (pGpC) base-pairs with a total of 715 
atoms. Periodic boundary conditions are used to simu- 
late infinitely long DNA chains well separated from each 
other, with negligible mutual interactions. Only the T 
/c-point of the 715-atom cell was included in the calcu- 
lations, which corresponds to 11 fc-points, if the helical 
symmetry is considered. The initial geometry was that 
of an experimental low-resolution A-helix structure for 
the pGpC sequence |3^, The geometry was relaxed by 
means of ab initio linear-scaling DFT, as described else- 
where The final structure had still well defined A- 
helix characteristics A second DNA realization 

(swapped hereafter) was built by swapping one single G- 
C pair into C-G out of the eleven in the cell, followed by 
another full relaxation. 

The electronic structure is analyzed in terms of the 
Kohn-Sham eigenvalues and eigenfunctions, based on 
Janak's theorem . They are obtained by a single cube- 
scaling diagonalization, using the final geometry and elec- 
tron density obtained with the linear-scaling method. 
The known band-gap problem of DFT is not important 
here, since the particular value of the band gap is of no 
relevance for the conclusions of this study as long as there 
is a gap, which is the case. The band widths have been 
shown to be well described by this method, at least up to 
the kind of accuracy needed in this work. A much more 
costly calculation including explicit electronic correlation 
would not be justified. Mott-like correlations, that could 
seem to be relevant for the narrow band-widths charac- 
teristic of the problem (see below), are not important 
given the small density of carriers (electrons or holes, 
see discussion below) expected in the system. The ef- 
fect of disorder and vibrations is explicitly considered 
below. The next section is devoted to a study of nucleic 
base pairs where the focus is mainly on the hydrogen 
bonds, which might in principle be worst described by 
our methodology. It will be followed by the results for 
A-DNA. 



A systematic study of 30 nucleic base pairs has been 
performed with the methodology presented above, in or- 
der to assess the reliability of the approximations, but 
also obtaining new results. For this study we have used 
a standard doble-C polarized (DZP) basis set, namely, a 
double (split valence) basis for each valence orbital plus 
polarization functions in all the atoms. The cutoff radii 
for the atomic orbitals of each element were obtained for 
an energy shift Q of 50 meV. 

A significant range of configurations of the four bases 
guanine, cytosine, adenine and thymine (G, C, A, T) are 
considered, the same as those studied by Sponer et al. 
in their MP2 study. The Watson-Crick configurations are 
designated WC, and the Hoogsteen, reversed Hoogsteen 
and reversed Watson- Crick appear as H, RH, RWC re- 
spectively. Other configurations are distinguished simply 
with numbers, eg. AAl, AA2, etc. following the nomen- 
clature of Hobza and Sandorfy Q . The structures of the 
bases and base-pairs studied in this work can be found 
in Figures 1 and 2 of Ref. |30| . For the numbering of 
the atoms we followed Ref. |44 In the following our re- 
sults are compared with those of other methods. The di- 
rect comparison with experimental geometries has been 
shown to be misleading because of the important charge 
rearrangements that happen in the crystal phases used 
in experiments p^ . 

The interaction energy Eint in the following is defined 
as the energy of the base-pair minus the energy of each 
base with the same geometry it has in the pair. The 
total stabilization energy Et is defined as the difference 
between the energy of the pair and that of each base in 
its isolated optimal geometry. The difference between 
both is thus the deformation energy, i.e., the increase in 
intramolecular energy due to the geometry change when 
the base-pair is formed. 

In this study all the energies have been corrected for 
basis set superposition error (BSSE) using the standard 
Boys-Bernardi counterpoise correction The correc- 
tion found for E^t is used for Et as well, given the fact 
that the relaxation of the isolated base molecules has 
been found to change the BSSE correction by less than 
10%. No BSSE correction was included in the forces. The 
effect of this approximation was gauged for the AAl base- 
pair, for which BSSE corrected and uncorrected forces 
were calculated for independent relaxations, showing a 
deviation of 0.02 A in the final acceptor-hydrogen dis- 
tances. 

Our results for Eint as compared with the MP2 re- 
sults of Sponer and coworkers pO| are presented in ta- 
ble I, as calculated for the same geometries obtained at 
the Hartree-Fock (HE) 6-31G** level @. For all the 
base-pairs except GG4, the agreement is quite satisfac- 
tory, with differences smaller than 8% and much less in 
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TABLE L Base- 


pair Interaction Energies 


( f^/A^+ in Wpal / TTiol 1 


at HF/6-31G** 


geometries. 






Pair 


Mr 2 


bi-tjo lA 


Deviation (%) 


GCWC 


-25.8 


-26.8 


-3.9 


GGl 


-24.7 


-25.1 


-1.6 


GCNEW 


-22.2 


-21.7 


2.2 


CC 


-18.8 


1 T C 

-17.5 


6.9 


GG3 


-17.8 


-16.6 


6.7 


GAl 


-15.2 


-15.5 


-2.0 


GTl 


-15.1 


-15.0 


0.7 


GT2 


-14.7 


-14.5 


1.4 


ACl 


-14.3 


-14.0 


2.1 


GCl 


-14.3 


-14.7 


-2.8 


AC2 


-14.1 


t A fr 

-14.7 


-4.2 


GAS 


-13.8 


-13.8 


0.0 


TAH 


-13.3 


-13.7 


-3.0 


TARH 


-13.2 


-13.6 


-3.0 


TAWC 


-12.4 


-12.3 


0.8 


TARWC 


-12.4 


-12.3 


0.8 


AAl 


-11.5 


1 1 ^7 

-11.7 


-1.7 


GA4 


-11.4 


-11.7 


-2.6 


TC2 


-11.6 


-10.8 


7.5 


TCI 


-11.4 


-10.6 


7.0 


AA2 


-11.0 


-11.4 


-3.6 


TT2 


-10.6 


-9.9 


6.6 


TTl 


-10.6 


-10.1 


4.7 


TT3 


-10.6 


-10.2 


3.8 


GA2 


-10.3 


-10.6 


-2.9 


GG4 


-10.0 


-7.4 


26.0 


AA3 


-9.8 


-9.8 


0.0 


2aminoAT 


-15.1 


-15.2 


-0.7 



"From ref. 30 



most cases. GG4 is an exception to the general trend 
with a relative difference of 26% that is reduced to 16% 
by extending the basis range to the longer radii given by 
an energy shift of 10 meV. A similar correction for the 
other base pairs gives much smaller variations. The stan- 
dard deviation of our results compared to the MP2 values 
is of 0.73 kcal/mol, the same value found by Sponer et 
aL|3^using the B3LYP functional |48| on the same geom- 
etry. The standard deviation between both sets of DFT 
results (our PBE and their B3LYP) is of 0.85 kcal/mol. 

Concerning geometries, the only pair relaxed with MP2 
is cytosine-cytosine [Q. Table ^ compares the length 
of the hydrogen bond and both Eint and Et for various 
methods, showing again a very satisfactory agreement. 

The results shown in this section validate the approxi- 
mations used in our method for the system under study, 
including the sensitive hydrogen bonds. Even if the pil- 
ing up of bases was not investigated, other studies show 



TABLE IL Cytosine-cytosine base-pair, distances in A and 
energies in kcal/mol. 

HF" MP2/HF'' MP2'' DFT-^ This work 



d(N4(H)-N3) 3.050 

Eint -17.3 
Et 



3.050 

-18.8 
-17.5 



2.980 2.900 



-20.5 
-18.7 



-20.4 
-18.1 



2.872 

-21.1 
-18.5 



"Energies and geometry with HF 6-31G** | po| 
''MP2 energies on HF geometry 
'^MP2 energies and geometry 
''B3LYP energies and geometry 



TABLE III: Hydrogen-bond distances (in A) for the guanine- 
cytosine pair. DZP and DZ(P) stand for two different basis 
sets used in this work. 







B3LYP" 


BP86'' 


VWN-BP" 


DZP 


DZ(P) 


N2 


- 02 


2.93 


2.87 


2.930 


2.872 


2.892 


N2 


- H 






1.035 


1.035 


1.023 


Nl 


- N3 


2.92 


2.88 


2.923 


2.913 


2.822 


Nl 


- H 






1.051 


1.056 


1.043 


06 


- N4 


2.78 


2.73 


2.785 


2.770 


2.715 


H - 


N4 






1.055 


1.057 


1.054 



"J. Sponer et al. ref. |30| 

''C. F. Guerra et al. ref. ^ 

"^R. Santamarfa and A. Vazquez, ref. |50| 



that GGA functionals can describe them with the re- 
quired accuracy when connected by the sugar-phosphate 
chains. 

The whole set of base pairs as well as the isolated 
bases have been relaxed with our method. The results 
for geometry and interaction energies can be found in 
pof . We extract in table III the results for the guanine- 
cytosine pair, of relevance for the DNA calculations in the 
next section. The results of several authors ^ |5^ 
are compared, which include different GGA function- 
al, namdy B3LYP @, BP86 [HI ||, and VWN-BP 



m ||. The results of the DZP basis of this study 
are also compared with those of the basis set used for 
DNA below, DZ(P), a double-C basis @, H], which was 
polarized for phosphorous and for the atoms involved in 
the hydrogen bonds. The Kohn-Sham eigenvalues around 
the Fermi level for an isolated GC pair were checked for 
DZ(P), which compare well with the corresponding val- 
ues for the DZP basis. The convergence of the spatial 
range of the basis orbitals, important for band widths in 
the DNA molecule, was also tested [|55|. 
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RESULTS FOR A-DNA 



Figure |] shows the bands of pGpC close to the Fermi 
level. As expected for a well saturated molecule, the 
Fermi level is in a clear band gap, which is of 2.0 eV in 
this case (this number is to be taken with caution since 
usual DFT functionals tend to underestimate the band 
gaps of insulators by a substantial amount). Electronic 
transport will thus be mediated by carriers in either the 
top of the valence band (holes, or radical cations in bio- 
chemical language) or in the bottom of the conduction 
band (electrons or radical anions), introduced by doping. 

In the absence of disorder, the electronic structure 
close to the Fermi level shows well defined bands with 
eleven states per unit cell, one per base pair in the unit 
cell. The top-most valence band has a very small band- 
width of 40 meV (see figure |l|), the energy separation 
with the next band below being ten times larger. This 
first band is associated to the 7r-like highest occupied 
molecular orbital (HOMO) of the guanines. The con- 
tinuous line in the figure is the tight-binding band of a 
one-dimensional system with one orbital per unit cell. 
The excellent agreement of this extremely simple model 
with the ab initio points [see figure |l| (b)] demonstrates 
that each base pair contributes to this band with one sin- 
gle orbital that interacts negligibly with other orbitals in 
the pair. This allows the modelling of the basic physics 
of transport in this molecule with one single orbital per 
base pair with an electron hopping interaction of 10 meV. 
Other terms for a model Hamiltonian are quantitatively 
justified below. 

Figure ^ displays the charge density of the states asso- 
ciated to this band, which appears almost exclusively on 
the guanines, with weight neither in the backbones nor in 
the cytosines. The lowest conduction band has a width of 
270 meV, similar to the width of the spectroscopic feature 
observed by Porath et al. It is separated from the 
next band by 0.7 eV. Similarly to the HOMO situation, 
this band is made of the lowest unoccupied molecular 
orbital (LUMO) of the cytosines, as shown in figure ||, 
although the overlap in this case is more pronounced, as 
corresponds to an eight times larger bandwidth. 

Note that any matrix element between HOMO and 
LUMO states is very much depressed by the very small 
spatial overlap between them. Optical absorption will 
thus be very weak for photon energies corresponding to 
the HOMO-LUMO gap (adequately corrected including 
correlation and excitonic effects). The unoccupied band 
above the LUMO band remains of a marked cytosine 
character. The third unoccupied band 3.9 eV above the 
HOMO band (see figure |l|) is the first one with substan- 
tial weight on the guanines, and thus with an apprecia- 
ble oscillator strength for absorption. The comparison of 
these results with those obtained for wet conditions ||2^ 
shows that the LUMO is quite sensitive to the environ- 
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FIG. 1: (a) Kohn-Sham energy bands (points) close to the 
Fermi level as functions of the k values corresponding to the 
helical symmetry. A indicates the HOMO-LUMO gap, Ai the 
intra-guanine gap, Wh the band width of the highest occupied 
band, and Wl the band width of the lowest unoccupied band, 
(b) Kohn-Sham highest occupied band (points), and band of 
a ID tight-binding model with one orbital per unit cell (line), 
with t = 10.1 meV, and t' = 1.5 meV as first and second 
nearest neighbor interactions, respectively. 
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FIG. 2: (Color stereogram) Surfaces of constant density (1.5 
10^ e~/A^) for the states corresponding to the lowest unoc- 
cupied band (red), and highest occupied band (blue) of the 
ordered pGpC structure. The graphs were produced with 
Molekel 



ment, moving from the cytosines to the cations when in 
presence of Na"*" and water. 

The effect of the swap on the electronic structure of 
the chain is dramatic, as can be seen in figure ^. The 
HOMO of the swapped guanine (indicated by an asterisk 
in figure 1) sinks 0.55 eV (14 times the HOMO band- 
with) into lower valence band levels. This stabilization 
is of electrostatic origin Figure ^ (left) shows the 

electronic density associated to the eleven highest occu- 
pied states for the swapped structure, showing the cut in 
the HOMO-state channel produced by the swapped pair. 
The situation is similar for the unoccupied band, albeit 
less dramatic, with the LUMO of the swapped cytosine 
raising 0.5 eV above the LUMO band. 

All the results presented so far refer to fixed, re- 
laxed geometries. It has been proposed, however, that 
the electron transfer could be mediated by phonons 
0, |l|, |l6[ |l|]. Indeed, there are many modes in DNA 
that are very soft and related to the nucleic bases, since 
these flat, rigid molecules hang from the backbone by a 
single bond. It is beyond the scope of this study to calcu- 
late phonon modes and their interactions with electrons, 
but interesting qualitative information on the matter can 
be inferred from the conjugate-gradient (CG) minimiza- 
tions. If close to the minimum, any snapshot during the 
minimization gives the atomic positions that correspond 
to the freezing of vibration modes with various ampli- 
tudes. The closer to the minimum, the more important 
the low-frequency modes versus the high-frequency ones 
p6[ . We study thus the electronic structure of geometries 
chosen from the CG path with an energy of around 0.9 
eV above the minimum, that would correspond to room 
temperature or lower p7|. The highest occupied energy 



-2.4 * 



> 



-2.8 - 



-3.2 



Regular 



Ragged 



Swapped 



* 



-n/a r n /a 



FIG. 3: Highest occupied electronic levels for the regular 
pGpC, the ragged, and the swapped structures. The aster- 
isk indicates a localized state with most of its weight on the 
HOMO of the swapped G. 



levels of such a geometry (called ragged) is presented in 
figure ||. If not as dramatic as the swapping, the effect 
of these "vibrations" on the electrons is also important, 
since the guanine HOMO one-particle states spread over 
an energy width eight times the original band width. 



DISCUSSION AND CONCLUSIONS 

Electronic transport in DNA requires the presence of 
carriers. In controlled biochemical experiments (in solu- 
tion), doping is induced by intercalating or attaching to 
the DNA specific molecules of known tendency to donate 
or accept electrons In the dry experiments intro- 

duced earlier there is uncontrolled doping of DNA. Since 
the negative charge of each phosphate group along both 
chains has to be compensated by a cation from solution, 
in the drying process cations precipitate along the DNA 
helix decorating it. In addition to a very probable disor- 
der in the cation positions and variability in the amount 
of water molecules hydrating them, there will also be 
deficiency or excess of counter-cations. Missing cations 
introduce holes in the valence band and cations in excess 
(or cations with higher charge) introduce electrons in the 
conduction band. The drying procedures of the dry ex- 
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FIG. 4: Isosurface as in figure for the top eleven occupied 
states of the swapped structure (left). The same for one par- 
ticular state (the 11th) for a density value 50 times smaller 
(right). 



periments do not allow estimating the concentration of 
carriers, not even their sign, even if (in view of the ion- 
ization potential of G and the electron affinity of C) the 
presence of holes has been argued to be more likely. It 
is reasonable to expect in any case that carrier concen- 
trations will remain low, and that for the more drastic 
drying procedures the likelihood is higher for the pres- 
ence of holes, since cations will be pulled away with the 
water. pGpC has been reported to be p-type in at least 
one experimental setup ]lO| . 

It has been shown above that the HOMO band width 
is one order of magnitude smaller than the energy scale 
associated to sequence disorder. In terms of a one- 
dimensional Anderson model, this leads to electronic lo- 
calization over very few base pairs. The localization is 
observed when plotting the separated eigenstates. Fig- 
ure^ (right) shows one electronic state within the HOMO 
energy range. The isosurface plotted is for a density 50 
times smaller that that on the left, showing clearly the 
localization of the state within a few base pairs. Se- 
quence disorder will therefore dramatically depress the 
possible band-structure-like DC conduction in DNA of 
arbitrary sequence (A-DNA). The possibility [p9[ of elec- 
tron hopping or tunneling, like in proteins |12[| is open, 
but it is not expected to show the conductances close to 
one quantum reported for A-DNA by some groups ^. 
For electron doping, disorder fluctuations are still com- 
parable or larger than the LUMO band width, and the 
localization of the electronic states is still considerable, 
albeit not as extreme as for the hole-doped case. 

It is known that certain types of disorder do not 
localize electrons as in the Anderson model. They rep- 
resent, however, very particular realizations that cannot 
be expected from A-DNA in general. Very interesting 



recent results [|6l| show possible electron delocalization 
in ID systems with more general long-range correlated 
disorder, as the disorder of natural DNA. More research 
is needed to see how robust is this result when includ- 
ing other effects like off-diagonal disorder or vibrations. 
In any case, the very large disorder-to-bandwidth ratio 
presented above makes extended states for holes very un- 
likely, the typical state being localized as displayed in 
figure ^ (right). 

The particular experiment of Porath et al. |9j used a 
controlled DNA sequence, namely pGpC, for which the 
sequence disorder arguments do not apply. They explain 
their results in terms of band-like conduction. The results 
and analysis of Berlin and coworkers propose band- 
like conduction in this particular system as well. It is 
tempting to support this band-like condution assuming it 
proceeds via electrons in the LUMO, given the agreement 
of the widths of our LUMO band and the feature they 
observe spectroscopically (the sensitivity of the LUMO 
to the environment mentioned above has to be kept in 
mind, however). Such an explanation would be by no 
means conclusive in the light of our "vibration" results, 
since a deterioration of conduction should be expected 
with temperature, which does not seem to happen. The 
spread of HOMO states with frozen-in vibrations indi- 
cates the energy scale of electron-phonon interactions, 
again very substantial as compared with the band widths. 
Even if a more rigorous study is needed to quantify these 
interactions and their effect, our rough estimations above 
allow us to expect a substantial reduction of conduction 
with temperature even below room temperature. 

Considering the floppy modes in DNA and their ef- 
fect on the electron states, and the fact that a hole on 
a guanine represents a positive charge amidst negatively 
charged guanines, the idea of self trapped (small) po- 
larons seems not at all unlikely. Polarons have already 
been suggested |l^ and supported [|o| in the litera- 
ture. Some of the models used have assumed the polarons 
to be of an origin similar to that found in polyacetylene, 
i.e. of off-diagonal (hopping) origin. Polarons of elec- 
trostatic origin are described by an electronic Hamilto- 
nian in which the diagonal matrix elements are the ones 
that change with inter-base separation and stabilize the 
polaron, as in Holstein's model [^. A purely classical 
electrostatic model based on the ab initio results for the 
geometry and ground-state electrostatic characteristics 
of the system displays polaronic charge excitations 
with a spread of a few base pairs. DFT calculations on 
the matter will be presented elsewhere p^ . 
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